#!/usr/bin/perl -w
use strict;
use lib "/home/x/Code/WEB_PHASTREGEX";
use DBI;
use UCSChooks;
# Jason Weirather 2010.
#this script will take any list of hits, and append 
#any feature you want to the right side of the list.
#it assumes that the list will start off with the
#chromosome start and finish of the hit.  After that, it will appened
#feature location and linking information.

# input MUST be ordered first by chromosome, then by start or it won't work.
# starting lines with '#' at their beginning are allowed and will be passed as comments

my $table_basename = 'encode';
my $feature_id = 107;
my $filtering = 'false';
my $max_distance = 0; # number of base pairs from feature to return the result
my $max_distance_from_midpoint = -1; # exclude big peaks.
if(scalar(@ARGV) < 2) { die "cat orderedoutputfile.txt | ./add_generic_feature.pl table_basename feature_id filtering_boolean max_distance_outside_feature max_distance_from_feature_midpoint\ntable_basename is the table from which you will draw this feature (ie feature_overlay or encode)\nfeature_id is the id of the feature in that database.\nfiltering is true or false.  if false, return all results reguardless of proximity to the feature, if true return only results near the feature.\n max_distance_outside_feature lets you return back results that are just near the feature but not necessisarilly inside of it.\nmax_distance_from_midpoint lets you specificy a maximum distance from the center of the peak that you will return results.  This is useful if you want to exclude very large peaks that may occur by mistake during liftover.\ntable_basename and feature_id are the only required fields.  by default filtering is false (off)\n"; }
$table_basename = shift @ARGV;
$feature_id = shift @ARGV;
if(scalar(@ARGV) > 0) { $filtering = shift @ARGV; }
if(scalar(@ARGV) > 0) { $max_distance = shift @ARGV; }
if(scalar(@ARGV) > 0) { $max_distance_from_midpoint = shift @ARGV; }

# output by print_matches.pl is an acceptable format.
# output by this program is also an acceeptable format

# lets start by bringing in the results from this file.
# we can do it one line at a time, and process it one chromosome at a time.

my $line;
my $previous_chromosome = '-1';
my @bstart;
my @bfinish;
my %chromhits;
$chromhits{'starts'} = \@bstart;
$chromhits{'finishes'} = \@bfinish;
my %chromosomes;
my $previous_start = -1;
my $chromosome = '-1';
my %appended;



while($line = <STDIN>) {
  ### pass any comments.
if($line=~/^#/) { print "$line"; }
else {
  chomp($line);
  my @fields = split(/\t/,$line);
  $chromosome = shift @fields;
  ## the other fields don't really matter right now but we will keep them all.
  if($chromosome ne $previous_chromosome && $previous_chromosome ne '-1') { # new chromosome!  
    if(defined($chromosomes{$chromosome})) { die "$chromosome chromosomes out of order. see notes on input order.\n"; }
    $chromosomes{$chromosome} = 1;
    ## push the buffer out here.  have to do it again after the while too.
      ####### read these in from 
    my ($feature_starts,$feature_finishes,$feature_peaks,$original_chromosomes,$original_starts,$original_finishes) = get_features($table_basename,$feature_id,$previous_chromosome);
    #for(my $i = 0; $i < scalar(@{$feature_starts}); $i++) { print "$chromosome\t$feature_starts->[$i]\t$feature_finishes->[$i]\n"; }
    process_phastcons($previous_chromosome,$chromhits{'starts'},$chromhits{'finishes'},$feature_starts,$feature_finishes,$feature_peaks,$original_chromosomes,$original_starts,$original_finishes,$filtering,\%appended,$feature_id);
    #do some initializing here.
    $previous_start = -1;
    my @blankstarts;
    my @blankfinishes;
    %appended = ();
    $chromhits{'starts'} = \@blankstarts; #initialize the buffer we are saving hits from this chromosome
    $chromhits{'finishes'} = \@blankfinishes; #initialize the buffer we are saving hits from this chromosome
  } 
  my $start = shift @fields;
  if($start < $previous_start) { die "starts are out of order. see notes on inputs (prev: $previous_chromosome $previous_start, this: $chromosome $start)!\n"; }
  else { $previous_start = $start; }
  my $finish = shift @fields;  
  my $ostring = '';
  foreach my $field (@fields) { $ostring .= "$field\t"; }
  chop($ostring); 
  $appended{"$chromosome\t$start\t$finish"} = $ostring;
  ### should have good checks for input
  push @{$chromhits{'starts'}}, $start;
  push @{$chromhits{'finishes'}},$finish;
  $previous_chromosome = $chromosome;
}
}

## flush the buffer one last time...
####### 
#my ($feature_starts,$feature_finishes,$feature_peaks,$original_starts,$original_finishes) = get_features($table_basename,$feature_id,$previous_chromosome);
my ($feature_starts,$feature_finishes,$feature_peaks,$original_chromosomes,$original_starts,$original_finishes) = get_features($table_basename,$feature_id,$previous_chromosome);
process_phastcons($previous_chromosome,$chromhits{'starts'},$chromhits{'finishes'},$feature_starts,$feature_finishes,$feature_peaks,$original_chromosomes,$original_starts,$original_finishes,$filtering,\%appended,$feature_id);
## chromosome should contain the most up to date chromosome still.

### at the end of the file lets add comments about what the feature we added is.
  my $username = 'x';
  my $tablebasename = '';
  my $zdatabase = "zebraphish";
  my $platform = "mysql";
  my $user = $username;
  my $zdbh = DBI->connect("dbi:$platform:$zdatabase",$user,'');
  my $zquery = "SELECT * FROM $table_basename" . '_info' . " WHERE db_id = $feature_id LIMIT 1";
  my $zst = $zdbh->prepare($zquery) or die "preparing query failed: $DBI::errstr";
  $zst->execute() or die "The execution of query field: $DBI::errstr";
  my $species;
  my $original_source;
  my $feature_name;
  my $cell_type;
  my $replicate;
  my $note;
  my $average_peak_height;
  my $average_peak_width;
  my $stddev_peak_height;
  my $stddev_peak_width;
  my $conversion;
  my $original_species;
  my $db_id;
  
  while(my $row = $zst->fetchrow_hashref()) {
    $species = $row->{'species'};
    $original_species = $row->{'original_species'};
    $feature_name = $row->{'feature_name'};
    $cell_type = $row->{'cell_type'};
    $replicate = $row->{'replicate_number'};
    $note = $row->{'note'};
    $replicate = $row->{'replicate_number'};
  }
  print "###Feature Boolean -> For feature: $feature_name, cell_type: $cell_type, replicate: $replicate, species: $species, original_species: $original_species, note; $note\n";
  print "###Feature Distance -> For feature: $feature_name, cell_type: $cell_type, replicate: $replicate, species: $species, original_species: $original_species, note; $note\n";
  print "###Feature Score -> For feature: $feature_name, cell_type: $cell_type, replicate: $replicate, species: $species, original_species: $original_species, note; $note\n";
  print "###Feature Location -> For feature: $feature_name, cell_type: $cell_type, replicate: $replicate, species: $species, original_species: $original_species, note; $note\n";
  print "###Feature Source Location -> For feature: $feature_name, cell_type: $cell_type, replicate: $replicate, species: $species, original_species: $original_species, note; $note\n";

### end main




sub get_features {
  my $table_basename = shift(@_);
  my $feature_id = shift(@_);
  my $chromosome = shift(@_);
  my @starts;
  my @finishes;
  my @scores;
  my @original_chromosomes;
  my @original_starts;
  my @original_finishes;
  my $username = 'x';
  my $tablebasename = '';
  my $zdatabase = "zebraphish";
  my $platform = "mysql";
  my $user = $username;
  my $zdbh = DBI->connect("dbi:$platform:$zdatabase",$user,'');
  my $zquery = "SELECT * FROM $table_basename WHERE chromosome = '$chromosome' AND feature_overlay_info_id = $feature_id ORDER BY chromosome, start";
  my $zst = $zdbh->prepare($zquery) or die "preparing query failed: $DBI::errstr";
  $zst->execute() or die "The execution of query field: $DBI::errstr";
  while(my $row = $zst->fetchrow_hashref()) {
    push @starts, $row->{'start'};
    push @finishes, $row->{'finish'};
    push @scores, $row->{'score'};
    push @original_chromosomes, $row->{'original_chromosome'};
    push @original_starts, $row->{'original_start'};
    push @original_finishes, $row->{'original_finish'};
  }
  return(\@starts,\@finishes,\@scores,\@original_chromosomes,\@original_starts,\@original_finishes);
}

sub process_phastcons {
  my $chromosome = shift @_;
  my @hit_starts = @{shift @_};
  my @hit_finishes = @{shift @_};
  my @starts = @{shift @_};
  my @finishes = @{shift @_};
  my @scores = @{shift @_};
  my @original_chromosomes = @{shift @_};
  my @original_starts = @{shift @_};
  my @original_finishes = @{shift @_};
  my $conservation_index = 0;
  my %conserved_ids;
  my $filtering = shift @_;
  my %appended = %{shift @_};
  my $feature_id = shift @_;
  #print "checking " . scalar(@hit_ids) . " hit ids\n";
  #print "checking against " . scalar(@starts) . " starts\n";
  if(scalar(@hit_starts) != scalar(@hit_finishes)) { die "something wrong\n"; }

  my $username = 'x';
  my $tablebasename = '';
  my $zdatabase = "zebraphish";
  my $platform = "mysql";
  my $user = $username;
  my $zdbh = DBI->connect("dbi:$platform:$zdatabase",$user,'');
  my $zquery = "SELECT * FROM $table_basename" . '_info' . " WHERE db_id = $feature_id LIMIT 1";
  my $zst = $zdbh->prepare($zquery) or die "preparing query failed: $DBI::errstr";
  $zst->execute() or die "The execution of query field: $DBI::errstr";
  my $species;
  my $original_species;
  while(my $row = $zst->fetchrow_hashref()) {
    $species = $row->{'species'};
    $original_species = $row->{'original_species'};
  }

  for(my $i = 0; $i < scalar(@hit_starts); $i++) {
    ### need to do all outputs in this loop to keep them in order.
    my $hit = 'false';
    my $j = $conservation_index;
    my $too_far = 'false';
    #print "$j\n";
    while($j < scalar(@starts) && $hit eq 'false') {
      if($hit_starts[$i] > $finishes[$j]) { $conservation_index = $j; }
      if($starts[$j] > $hit_finishes[$i]) { $hit = 'skip'; }
      #print "trying $starts[$j] $finishes[$j] $hit_starts[$i] $hit_finishes[$i]\n";
      if($hit_starts[$i] >= $starts[$j] && $hit_finishes[$i] <= $finishes[$j]) {
        #print "true\n";
        $hit = 'true';
      }
      if($hit_starts[$i] <= $starts[$j] && $hit_finishes[$i] >= $starts[$j]) {
        #print "true\n";
        $hit = 'true';
      }
      if($hit_starts[$i] <= $finishes[$j] && $hit_finishes[$i] >= $finishes[$j]) {
        #print "true\n";
        $hit = 'true';
      }
      if($hit_starts[$i] <= $starts[$j] && $hit_finishes[$i] >= $finishes[$j]) {
        #print "true\n";
        $hit = 'true';
      }
      if($hit eq 'true') { 
        #### this is where it is a result.
        my $feature_midpoint = int(($finishes[$j]-$starts[$j])/2)+$starts[$j];
        my $hit_midpoint = int(($hit_finishes[$i]-$hit_starts[$i])/2)+$hit_starts[$i];
        my $distance = abs($feature_midpoint-$hit_midpoint);
        if($max_distance_from_midpoint == -1 || $distance < $max_distance_from_midpoint) {
          $conserved_ids{$i} = 1; # do the other ones later 
          ### print out a conserved result
          my $ostring = $appended{"$chromosome\t$hit_starts[$i]\t$hit_finishes[$i]"};
          print "$chromosome\t$hit_starts[$i]\t$hit_finishes[$i]\t$ostring\t";
          my $feature_url = UCSChooks->make_excel_link($species,$chromosome,$starts[$j],$finishes[$j]);
          my $original_url = UCSChooks->make_excel_link($species,$original_chromosomes[$j],$original_starts[$j],$original_finishes[$j]);
          print "1\t$distance\t$scores[$j]\t$feature_url\t$original_url\n";  ## add 4 columns
        } else { $hit = 'false'; }
      }
      $j++;
    }
    if($hit ne 'true' && $filtering ne 'true') { # we need to output here if
      my $ostring = $appended{"$chromosome\t$hit_starts[$i]\t$hit_finishes[$i]"};
      print "$chromosome\t$hit_starts[$i]\t$hit_finishes[$i]\t$ostring\t";
      print "\t\t\t\t\n";
    }
  }
  #print scalar(@conserved_ids) . " conserved ids\n";
  ### Each index of conserved_ids has the index of a conserved hit.
  # if we are filtering, don't annotate everything, just the ones with hits
  #if($filtering ne 'true') {
  #  for(my $i=0; $i <  scalar(@hit_starts); $i++) {
  #    if(!exists($conserved_ids{$i})) {
  #      # that should be the original.
  #    }
  #  }
  #}
}

